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Abstract 

We present the first two-dimensional general relativistic (GR) simulations of stellar core collapse 
and explosion with the CoCoNuT hydrodynamics code in combination with the Vertex solver for 
energy-dependent, three-flavor neutrino transport, using the extended conformal flatness condition for 
approximating the spacetime metric and a ray-by-ray-plus ansatz to tackle the multi-dimensionality of 
the transport. For both of the investigated 11.2 and ISA/q progenitors we obtain successful, though 
seemingly marginal, neutrino-driven supernova explosions. This outcome and the time evolution of the 
models basically agree with results previously obtained with the Prometheus hydro solver includ- 
ing an approximative treatment of relativistic effects by a modified Newtonian potential. However, 
GR models exhibit subtle differences in the neutrinospheric conditions compared to Newtonian and 
pseudo-Newtonian simulations. These differences lead to significantly higher luminosities and mean 
energies of the radiated electron neutrinos and antineutrinos and therefore to larger energy-deposition 
rates and heating efficiencies in the gain layer with favorable consequences for strong nonradial mass 
motions and ultimately for an explosion. Moreover, energy transfer to the stellar medium around the 
neutrinospheres through nucleon recoil in scattering reactions of heavy-lepton neutrinos also enhances 
the mentioned effects. Together with previous pseudo-Newtonian models the presented relativistic 
calculations suggest that the treatment of gravity and energy-exchanging neutrino interactions can 
make differences of even 50-100% in some quantities and is likely to contribute to a finally successful 
explosion mechanism on no minor level than hydrodynamical differences between different dimensions. 

Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics — relativity 
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1. INTRODUCTION 

More than 40 years after the first attempts by 
iColgate fc White! (|1966( 1. the quest for the supernova 
explosion mechanism is still ongoing. Since the solu- 
tion now seems (once again) within reach as several 
groups hav e come up with successful explosion models 
(iMarek fc Janka 2009; S.u _ wa et al. 2010; T akiwa ki et al.l 
1201 It iBruenn et aLf^OOgl : iBurrows et al.llmel l2007aD . 
the demand for accurate predictions of the neutrino and 
gravitational wave signals and of the nucleosynthetic 
yields becomes more urgent and naturally requires in- 
creased levels of sophistication in numerical core-collapse 
supernova simulations. 

Currently, there is a broad consensus that multi- 
dimensional hydrodynamical instabilities are one of 
the pivotal elements of the supernova problem. It 
has long been recognized that convection in the 
hot-bubble region serves to increase t he efficiency 
of neutrino h eating behind th e shock (iHerant et all 
T992. 1994; B urrows et al.l [19951 : I.Tanka fc MiilleH 119961 : 
MiiUer fc Jankal Il997n . and that another more re- 
cently discovered large-scale instabilit y, the stand- 
ing accretion-shock instability ("SASI", Blondin et alj 
2003; Blondin & Mczzacappa 2006; Foghzzo et al. 200l 
Ohnishi et al. 2006; Foglizzo et al. 2007; Schcc k et ail 
"OO&Iwakami et al...200&..2009. ; Eernandez k Thompson! 
20091 : iFbrnandezI 120101 ) has a similar beneficial effect. 



Both instabilities help to keep the accreted material in 
the gain region for a longer time before it is advected 
deeper into the cooling region and onto the neutron star 



surface (|Buras et al.l l2006al: iMurphv fc Burrows! [20081) . 
If the advection time-scale Tadv through the gain region 
(sometimes also termed "residence time" of the matter 
in the gain region) is increased sufficiently and becomes 
comparable to or larger than the heating time-scale Thcat 
required to unbind the material between gain radius and 
shock, a runaway situation occurs, in which neutrino 
heating leads to shock expansion, which in turn length- 
ens the residence ti me, thus again increasing t he neutrino 
heating efficiency ( Janka fc K cil 199S; iThomps on 20001 
Janka et al.! !2001l; Buras et al.. .2006al : !Thompson et all 
2005; Mur phy fc Burrows!!2008l ) . 

However, whether an explosion can actually be brought 
underway in this fashion in the most sophisticated su- 
pernova models with detailed neutrino transport has 
not yet been finally and unambiguously established on 
the basis of state-of-the-art neutrino hydrodynamics 
simulations in 2D. Using their Vertex-Prometheus 
code, which employs a variable Eddington factor tech- 
nique to solve the neutrino moment equations and 
the "ray-by-ray-plus" approach to cope with multi- 
dimensional transport, the Garching group has found 
explosions by the SASI-aided neutrino-dr iven mecha- 
nism for an il.2M(^ solar rn ass progenitor (jBuras et al.l 
!2006at !Marek fc Jankal [20091) . which can be reproduced 
robustly for stiffer and softer nuclear equations of state, 
(Marek et al. in preparation) and for a 15AfQ progenitor 
with artificially imposed rotation at a rather late time 
w 550 ms after bounce. In other cases, i.e. for the same 
15AfQ progenitor without rotation and with a better ef- 
fective gravitational potential, and for more massive pro- 



2 



Miiller, Janka & Marek 



genitors, no explosion was observed until 400 — 500 ms 
after bounce. In contrast to this, the Oak Ridge group, 
relying on a multi-group flux-limited diffusion (MGFLD) 
algorithm combined with the ray-by-ray-plus approach, 
has obtained explosions for a host of different progeni- 
tors (Bruenn et al. 20061 I2009D . while the former Arizona 
group did not obtain neutrino-driven explosions at all 
with their 2D MGFLD code Vulcan (jLivne et al.ll2004 
|2007[) but found acoustic explosions powered by sound 
waves create d by vi olent dipolar oscillations of the proto- 
neutron star ([Burro ws et al. 2006, 2007b). Using yet an- 
other neutrino transpor t scheme, the "isotropic d iffusion 
source approximation" ([Liebendorfer et aT]|2009[ ) in the 
ray-by-ray-approximation and without a nd Vr trans- 
port, an ex plosion has also been reported by iSuwa et al.l 
(poTol ) and iTakiwaki et al.l (poTll ) for an II.2M0 and a 
I3M0 progenitor. 

Given the disparity of methods and input physics, 
these different results should not be overly surprising: 
None of the aforementioned groups follow approaches 
that are completely identical with respect to the treat- 
ment of neutrino transport (variable Eddington factor 
method vs. MGFLD vs. IDS A, inclusion/omission of 
non-conservative scattering, velocity effects, and gravi- 
tational redshift, ray-by-ray transport vs. multi-angle 
transport), the neutrino processe s (e. g. omission of 
u/r neutrinos in iSuwa et all 120101 and ITakiwaki et al.l 
1201 1[) . the hydro solver (high- resolution shock captur- 
ing schemes vs. artificial viscosity), or the treatment of 
gravity (effective relativistic potential vs. Newtonian ap- 
proximation) and nuclear burning (network vs. flashing 
vs. omission of burning) . A clear sensitivity of supernova 
simulations to Newtonian vs. GR gravity, the sophistica- 
tion of the neutrino opacities, and observer corrections in 
the transport equ ations has re c ently been demonstrated 
in ID models by iLentz et al.l (|2011[ ). Different results 
depending on the input physics and the approximations 
used by the different groups are all the more to be ex- 
pected considering that some of the 2D models, whether 
exploding or non-exploding, appear to be marginal cases 
(see e.g. IMarek &: J anka 2009) anyway. From this per- 
spective, the viability of the neutrino-driven mechanism 
in 2D remains a controversial issue and should still be 
investigated further with the help of better simulations. 

Considering the status of current 2D core-collapse 
supernova models, it is conceivable that they still 
miss a crucial ingredient for robust explosions. Al- 
ternatives to the standard neutrino-driven mecha- 
nism have therefore been proposed and explored, 
such as magnetohydro dynamically-driven explosions (see 
iBurrows et al.|[2b07 a' and references therein), energy in- 
put by acoustic waves (Burrows et al. 2006, 2007b), or 
a phase transition to quark matter (jSagert et al.l 120091 ). 
Moreover, 3D effec ts have recently b een ad vocated as the 
decisive factor by iNordhaus et al.l ()2010[ ) in their com- 
parison of ID, 2D, and 3D simulations using a strongly 
simplified prescription for neutrino heating and cool- 
ing. Adopting the conceptual view of IBurrows &: Goshvl 
(jT993) of the neutrino-driven mechanism as a critical 
phenomenon, they report a reduction of the so-called 
critical luminosity in 3D by 15% - 25% compared to 2D, 
and single out the dimensionality as the key to successful 
supernova explosions, discounting other factors such as 
general relativity and detailed neutrino microphysics as 



minuscule corrections. However, whether these results, 
obtained by means of a very rough approximation for 
neutrino heating and cooling, can be verified by simula- 
tion s with an elaborate transport treatment i s yet unclear 
(see lHanke et al.l l2011: Takiwaki et al. 2011). Moreover, 
the lower estimate of Nordhaus et al, (2010,) for the im- 
portance of 3D effects relative to 2D would probably not 
make them the single most important factor in supernova 
physics, at least not by far. While there is no doubt that 
3D models are indispensable for better understanding 
the S ASI (Blondin & Mezza cappa 2007; Iwakami et al] 
12001 12009t lFernlndeg[20lOD . the morphology of super- 



nova explosions, and the kicks an d spins of neutron 
stars ([Hammer et al.ll20ToHW ongwa thanarat et al.ll2010l : 
iRantsiou et al.ll2011l: lFernande^'l2010[ r our understand- 
ing of supernovae certainly does not hinge on dimension- 
ality alone. 

General relativity (GR) is undoubtedly another ma- 
jor factor in supernova physics due to the compact- 
ness of the neutron star {M/R sa 0.1... 0.2 in ge- 
ometrized units) and the occurrence of large infall (w/c « 
—0.15... — 0.3) and outfiow velocities. In spherical 
symmetry (ID), the gauge freedom can be exploited to 
take GR effects into account without sacrificing accu- 
racy in the neutrino transport sector more readily than 
in 2D and 3D; gray a nd multi-group flux-limited dif- 
fusion schemes in ID ([Baron et al.l Il989t iBruenn et al.l 
120011 ) have therefore been available since the 1980s and 
were followed b y general relati vistic Boltzmann solvers 
a few years ago ([Yamada et al.llT999,: .Liebendorfer et al.l 
I2001ll2004[). Comparis ons with the Newtonian case in 
ID (I Bruenn et aO2001f) clearly showed the importance of 
GR effects, p articularly at late tim es during the accretion 
phase, where iBruenn et al.l ([20011 ) find a much stronger 
recession of the shock (50% smaller radius), higher lumi- 
nosities, in particular for electron neutrinos and antineu- 
trinos (by ~ 40%) and higher mean neutrino energies (by 
sa 15%) in the GR case. This strong sensitivity to GR ef- 
fects was recently c onfirmed with up-to-date neutrino mi- 
crophysics by Lent z et al.l ([20111 ). Considering the mag- 
nitude of relativistic effects, it is obvious that they ought 
to be properly considered in multi-dimensional supernova 
models as well. 

Unfortunately, until recently the only viable op- 
tion to include GR effects in multi-dimensional neu- 
trino hydrodynamics simulations while retaining the 
framework of Newtonian hydrodynamics has been the 
use of modified gravitational or "effective" potentials 
based on the Tol man-Oppenheimer - Volko v equation of 
stellar structure ( Rampp fc Jankal 120021 : IMarek et al.l 
120061 iMiiller et al.l boOSD. This "pseudo-Newtonian" ap- 
proach is ea sy to implement, yields very satisfactory re- 
sults in ID CLiebe ndorfer et all 12005: Ma rek et al.l 120061: 
IMiiller et al..,2010i) . and certainly provides a rough first- 
order approximation for GR effects in multi-dimensional 
supernova models. As a complementary approach start- 
ing from multi-dimensional general relativistic hydro- 
dynamics simulations, there have been efforts to ad- 
dress certain aspects of core-collapse supernovae with the 
help of heavily simplified neutrino treatment s such as a 
"deleptonization scheme" ([Liebendorfeiil2005h or a trap- 
ping scheme, but the applicability of such methods is 
naturally limi ted, e.g. to the colla pse and bounce phase, 
(|Ott et al.. i2007at iDimmelmeieT 3t al.. ,2007a ,2008,), or 
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fast black hole formation (|Ott at al.ll2011[ ). 

In this paper, we pursue a considerably more ambi- 
tious course and present the first 2D core-collapse super- 
nova models combining general relativistic hydrodynam- 
ics and a sophisticated energy- dependent neutrino trans- 
port scheme. Using the Vertex-CoCoNuT code with 
a ray-by-ray-plus treatment of multi-flavor, multi-energy 
2D ne utrino transport as documented in iMiiller et alj 
(|2010l henceforth Paper I), we have conducted simula- 
tions of two progenitor models with II.2M0 and 15Mq 
beyond the onset of the explosion. These simulations are 
complemented by three non-exploding runs of the 15Mq 
progenitor, viz. one with a pseudo-Newtonian (effective 
potential) treatment of gravity, one with purely Newto- 
nian gravity, and one with a simplified set of neutrino 
interaction rates. The major purpose of these simula- 
tions can be summarized as follows, 

1. to demonstrate the feasibility of full-scale multi- 
dimensional general relativistic supernova simula- 
tions, 

2. to evaluate the role of general relativistic effects 
in the explosion mechanism and the quality of the 
pseudo-Newtonian approach, 

3. to gauge the sensitivity of the heating conditions 
to the neutrino physics input, 

4. to closely reproduce and verify pseudo-Newtonian 
simulations with the Prometheus-Vert ex code 
(|Buras et al.l [20063 IMarek fc Jankall2009D using a 
different hydrodynamics solver, thus eliminating 
reservations about the reliability of the numerics 
in existing 2D supernova models. 

Naturally, the impact of GR on the gravitational wave 
(and also on the neutrino) signal from core-collapse su- 
pernovae is also a question of paramount importance. 
However, this issue cannot be discussed within the scope 
of the present paper with its focus on the dynamics of 
the post-bounce evolution, and will be the subject of a 
follow-up publication. 

To address the aforementioned issues, we have struc- 
tured our paper in the following manner: We review the 
numerical scheme and outline the model setup and in- 
put physics used in this work in Sections [2] and [3l In 
Section HI we describe the shock evolution and the ex- 
plosion morphology of our relativistic explosion models. 
In Section [SJ we then address differences to the non- 
exploding runs and provide a more quantitative analysis 
of the heating conditions for our models in order to deter- 
mine the impact of the GR treatment and the neutrino 
microphysics on the post-bounce dynamics. We summa- 
rize the main results from this analysis in Section |6] and 
evaluate the implications with respect to the major issues 
of our paper. 

2. NUMERICAL METHODS 

We perform numerical simulations with the gen- 
eral relativistic neutrino hydrodyna mics code Vertex - 
CoCoNuT introduced in Paper I (jMiiUer et al.|[20Tol ). 
which is a combination of the neu t rino transpo r t solver 
Vertex (Rampp i^Jankal 120021: iBuras et al.l l2006bD 
and the relativistic hydrodynamics code CoCoNuT 



(|Dimmelmeier et al.|[2002ll2005D . CoCoNuT is a time- 
explicit, directionally unsplit Eulerian Godunov-type 
finite-volume solver written for spherical polar coordi- 
nates and uses piecewise parabolic (PPM) reconstruc- 
tion and Runge-Kutta time-stepping to achieve higher- 
order spatial and temporal accuracy. Our implementa- 
tion of CoCoNuT employs a relativi stic version of the 
HLLC approximative Riemann solver (jMignone &: Bodol 
|2005[) , but ad aptively switc hes to the more diffusive 
HLLE solver ([Einfeldtl I1988D in the vicinity of strong 
shocks to av oid the phenomenon of odd-even-decoupling 
(|Quirklll994D . In order to reduce spurious numerical en- 
ergy generation, we use an improved formulation of the 
energy equation described in Paper I. The metric equa- 
tions are solved approximatively using the extended con- 
forma l flatness condition (xCFC, ICordero-Carrion et al.l 
|2009() . a constrained scheme with improved numerical 
stability properties co mpared t o the original conformally 
flat approximation of llsenber^ ([1978). In the context of 
core-collapse, this approximation works extremely well as 
demonstrated by the excellent agreement between rota- 
tional core-collapse simulations using the CFC approx- 
imation and thos e relying on the full ADM formalism 
(|Ott et al.ll2007allbl : IDimmelmeier et al.|[2007bl) . 

The time-implicit neutrino transport module Vertex 
integrates the moment equations for the neutrino energy 
and momentu m density usin g a variable Eddington fac- 
tor technique (IRampp fc Jan ka 2002). The higher mo- 
ments of the neutrino radiation field that are required to 
close the system of moment equations are obtained from 
a simplified Boltzmann equation that is solved in con- 
junction with the neutrino moment equations within a 
fixed-point iteration. All velocity- and metric-dependent 
terms are fully included in the moment equations, as 
is energy redistribution by non-conservative scattering. 
In 2D, we resort to t he "ray-by-ray-plus " approximation 
(jBuras et al.ll2006bl : iBruenn et all 120061) . assuming that 
the neutrino distribution function is axially symmetric 
around the radial direction (which implies a radial fiux 
vector). This allows us to solve independent ID trans- 
port problems along "rays" corresponding to the angular 
zones of the polar grid. However, the lateral advection 
of neutrinos and the effect of lateral neutrino pressure 
gradients are both included to avoid unphysical b ehavior 
in the optically thick regime (jBuras "etani2006bD . 

We use an up-to-date set of neutrino interaction rates 
in Vertex; in addi tion to the "standard" set of neutrino 
opacities from iBruennI (Il985f). we take nucl eon-nucleon 
bremsstrahlung ([Hannestad fc RaffellJ I1998D . neutrino- 
neutrino pair conversion ( Buras et al.l l2003f) and inelas- 



tic ne utrino scattering off heavy nuclei (|Langanke et al.l 
I2OO8D into account. Furthermore, we compute elec- 
tron capt ures on heavy n u clei us ing the improved rate 
tables of iLanganke etall (|2003| ) instead of the FEN 
rates (iFuUer et al.lll982D . and as an alternative to the 
"conservative" (i.e. no energy exchange with nucleons 
treate d as vacuum particles) approximation of IBruennI 
1985f) we inclu de recoil effects, high-density correlations 
Burrows fc Sawve r 1 998, 1999), and weak magnetism 



corrections ( Horowit^ 1 1 99% in our treatment of charged 
and neutral-current neutrino-nucleon interactions (see 
Table 121). 
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3. MODEL SETUP 
3.1. Progenitors and Neutrino Physics 

In this paper, we consider two differei it non-rotating 
progenitors, namely t he models si 1.2 of iWooslev et aLl 
|Ib02) and sl5s7b2 of IWooslev fc WeaveI^ (|1995n . These 
11.2M0 and ISM© progenitors have been chosen such 
as to facilitate a compariso n with t he ps eudo-Newtonian 
simula tions of iBuras et al.l (|2006aD and iMarek fc Jan^ 
(120091). The ISM© ca se is of particular interest since 
IMarek fc Jankal ()2009[ ) observed an explosion for this 
progenitor in a run with artificially imposed rotation and 
an overly strong effective gravitational potential (their 
model LS-rot) while the corresponding non-rotating 
model (L S-2D), computed w ith the best effective po- 
tential of IMarek et al.l ()2006D . failed to explode before 
the end of their high-resolution run 420 ms aft er bounce. 
Within the model assumptions of M arek fc Jan ka (2009), 
the progenitor sl5s7b2 is therefore a marginal case, and 
hence ideally suited for investigating possible effects of 
general relativity and testing the accuracy of the effec- 
tive potential approach. In addition to two relativistic 
models (Gil and G15) for the two progenitors computed 
with Vertex-CoCoNuT, we therefore also consider a 
purely Newtonian (N15) and a. pseudo-Newtonian (M15) 
simulation of the ISM© star of lWooslev fc Weaver] (jl995| ) 
computed with Vertex-Prometheus. Furthermore, 
we also include another relativistic calculation (SI 5) with 
a slightly simplified set of neutrino opacities to assess the 
importance of improved interactions rates (particularly 
for neutrino- nucleon processes) . A summary of the input 
physics for these models is given in Table [l] and the dif- 
ferences between the full and simplified set of neutrino 
rates are given in detail in Table [5J 

With the exception of model N15 (where 128 angular 
zones were used), all runs were performed on a spherical 
polar grid covering 180° in latitude with 400 x 64 zones 
initially. The distribution of the radial grid varies be- 
tween the individual simulation^, but was chosen (and 
when necessary re-adjusted) such that the density gra- 
dient at the surface of the proto-neutron star remains 
sufficiently well resolved during the simulation. With a 
moderate angular resolution of 2.8° we have settled for 
an affordable compromise in the first long-time multi- 
dimensional simulations using the relativistic Vertex- 
CoCoNuT code. For future calculations (relying on 
more efficient parallelization of the hydrodynamics and 
metric solvers), a higher resolution is clearly desirable, 
but past experience with pseudo-Newtonian simulations 
at least suggests that 64 angular zones already provide 
a reasonable level of accuracy to study systematic differ- 
ences between models with different input physics. Both 
for the 11.2Mq and the ISM© progenitor, this resolution 
is in fact sufficient for good quantitative agreement with 
high-resolution studies in the pseudo-Newtonian case as 
demonstrated b y the similarity of our model M15 with 
model LS-2D of IMarek fc 1^512001 For comparisons 
or extrapolation, it should nonetheless be borne in mind 
that runs with higher angular resolutio n generally appear 
somewhat more optimistic in 2D (cp. iHanke et al.ll201ll 

^ Specifically, stability considerations require different zoning in 
the hydrodynamics modules CoCoNuT and Prometheus near the 
origin of the spherical polar grid. 



for numerical tests with simple neutrino source terms). 
For the neutrino transport, we chose a logarithmically- 
spaced energy grid with 12 bins ranging from MeV to 
380 MeV, which - though not an optimal resolution - 
yields fairly similar dynamics compared to a finer zoning 
in energy space (cp. iMarek fc Janka.200& for an example 
with only 9 bins). 

3.2. Equation of state 

The simulations were perfor med using a soft version o f 
the equation of state (EoS) of iLattimer fc Swestvl (|1991l ) 
with a bulk incompressibility modulus of nuclear mat- 
ter of if = 180 MeV (LS180). The value of K in this 
EoS has raised some objections ()Nordhaus et al.l 120101 : 
iDasgupta et al.l 120 lit ) because measurements point to 
K = 240 MeV for symmetric nuclear matter, and be- 
cause observations have recently pushe d the maximum 
neutr on star mass to at least 2M© (jPemorest et al.l 
|2010[) . which is incompatible with this particular EoS 
(see Figure [T|) . 

However, this criticism neglects that (at least for the 
range of progenitors we are considering here) the crucial 
quantity for the development of the explosion is not the 
maximum neutron star mass, but rather the radius of 
the proto-neutron star during the reheating phase. For 
the baryonic neutron star masses we obtain at the end 
of our longest-running simulations (1.35M© and 1.58M© 
for the 11.2AfQ and 15A/© progenitor, respectively), the 
LS180 EoS yields a final radius (in the cold state) of 
around 12 km (see Figure [IJ, which is smaller by only 
0.5 ... 0.7 km than the neutron star radius predicted by 
the ("stiffer") LS220 EoS (which has K = 220 MeV and 
yields a maximum neutron star mass larger than 2Mq). 
Figure [T] also shows that the deviation is of similar mag- 
nitude for hot proto-neutron stars. Moreover, a neutron 
star radius of 12 km is well compatible with observed 
neutron star radii; it actually agrees nicely with the best 
existing (observational and theoretical) e stimates of neu- 
tron star radii, see e .g. th e papers bv Stei ner et al.l ()2010l ) 
and lHebeler et all (120101), a fact th at, for example, dis- 
favors the EoS of iShen et al.l (|1998| ). because the latter 
yields a radius of ~ 15 km for a neutron star with a grav- 
itational mass of 1.4M©. The prospective differences to 
a simulation with LS220 become even smaller if we con- 
sider that the two neutron stars are even less massive 
when the explosion develops (e.g. at a baryonic mass of 
1.51M© for the 15M© progenitor). 

It is also important to note that the neutron stars 
formed in the 11. 2M© and 15Mq runs have baryonic 
masses of 1.35M© and 1.58Mq, but the corresponding 
gravitational masses are about 10% lower than these 
values. These masses are far below the maximum gravi- 
tational mass that can be supported by the LS180 EOS 
(which is about 1.85Mq, corresponding to a baryonic 
mass well beyond 2Mq; Figure [IJ. 

In the neutron star mass regime we are dealing with in 
this paper, LS180 therefore remains a justifiable choice 
for the EoS, which yields neutron star radii quite similar 
to the "stiffer" version LS220, and is also very similar 
in most other EoS properties (see, e.g., the pressure in 
Figured]); sizable differences appear only close to the lim- 
iting neutron star masses supported by these equations 
of state. 
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Table 1 

Model Setup 







neutrino 


hydro 


treatment of 


final post-bounce 


angular 


model 


progenitor 


opacities 


solver 


relativity 


time reached 


resolution 


Gil 


sll.2 


full set 


CoCoNuT 


OR hydro + xCFC 


920 ms 


2.8° 


G15 


sl5s7b2 


full set 


CoCoNuT 


OR hydro + xCFC 


775 ms 


2.8° 


S15 


sl5s7b2 


reduced set 


CoCoNuT 


OR hydro + xCFC 


474 ms 


2.8° 


M15 


sl5s7b2 


full set 


Prometheus 


Newtonian + modified potential^ 


517 ms 


2.8° 


N15 


sl5s7b2 


full set 


Prometheus 


Newtonian (purely) 


525 ms 


1.4° 



''Case A of lMarek et"all 1200^ 1 



Table 2 

Neutrino physics input 





full rates 


simplified rates 


process 


(Gil, G15, M15, N15) 


(S15) 


uA ^ uA 


Horowitz (1997) Hon-ion correlations) 


Bruenn & Mezzacappa (1997) 




Langanke et al. (2008) (inelastic contribution) 


Itoh et al. (2004) (ion-ion-correlations) 


u ^ u 


Mezzacappa & Bruenn (1993) 


Mezzacappa & Bruenn (1993) 


u N N 


Burrows & Sawyer (1998)^ 


Bruenn (1985) 


Ve n ^ e~ p 


Burrows & Sawyer (1998)=' 


Bruenn (1985) 


Pe P ^ n 


mirrows & Sawyer ( 1998)=^ 


a-uenn (1985) 


Ue A' ^ e- A 


Langanke ct al. (2003) 


Fuller et al. (1982); Bruenn (1985) 


vv ^ e~ e+ 


Bruenn (1985); Pons et al"(1998) 


Bruenn (1985); Pons ct al. (1998J 


vv NN ^ NN 


Hannestad k. Raffelt (1998) 


Hannestad & Raffelt ( 1998) 




Buras et al. (2003) 




(-) (-) _{-) (-) 


Buraa et al. (2003) 





" Note that these reaction rates account for nucleon thermal motions, phase-space blocking, energy transfer to the nucleon associated with 
recoil ("non-conservative" or "non-isoenergetic " scattering), and nucleo n correlations at high densities. Moreover, we include the qu enching 
of the axial- vector coupl ing at h igh densities l lCarter Prakashll2003) . correction to the effective nucleon mass IIReddv et al.lfT999 ). and 
weak magnetism effects IjHorowitz 2002 ) in our full set of rates. 




R [km] R [km] 

Figure 1. Mass-radius relations of the equations of state LS180 (blue) and LS220 (red) for the gravitational mass (left panel) and the 
baryonic mass (right panel). Solid lines display the case of cold neutron stars (T = 0), while curves for the case of a hot proto-neutron star 
with a constant entropy of s = 1.5 k;, /nucleon are shown as dashed lines. The black horizontal line in the left panel corresponds to a mass 
of 1.97M0 as measured by IDemorest et al.l (|201(]l ) for the pulsar J1614-2230. The gravitational masses for neutron stars with baryonic 
masses of 1.36Mq and I.SSMq are indicated both for T = (solid blue horizontal lines) or s = 1.5 ki,/nucleon (dashed blue horizontal 
lines) in the left panel (figures provided by A. Bauswein.) 
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Figure 2. Pressure as a function of density for LS180 (blue) 
and LS220 (red) for four different temperatures (1 MeV, 5 MeV, 
10 MeV, 15 MeV from bottom to top) for Ye as in cold beta- 
equilibrium (figure provided by A. Bauswein). 




time after bounce [s] 



Figure 3. Average shock radius and proto-neutron star (PNS) 
radius (defined by a fiducial density of 10^^ g cm~'^) for the 2E) 
models G15 (GR, full rates, black thick solid line), S15 (GR, re- 
duced rates, blue, thick, dash-dotted), M15 (pseudo-Newtonian, 
full rates, red, thick, dashed), and M15 (purely Newtonian, black, 
thick, dotted). ID models corresponding to G15, M15, and S15 
are also shown as thin lines for comparison. Note that the shock 
is located considerably further out in S15-1D than in G15-1D and 
M15-1D. This is a consequence of the strong sensitivity of the shock 
position rah to the PNS radius, rpNs, for a stationary spherical ac- 
cretion flow (rsh DC rpVjo, see, e.g.. Equation (1) of Marek & Janka 
l2009iV The larger PNS radius in S15-1D can in turn be traced 
to less efficient cooling by /i/r neutrinos and higher tempera- 
tures in the density region 

10^2 , lo^^^ g cm ^. Different PNS 
radii (caused by proto-neutron star convection; see Appendix C in 
[Bur as et al. 2006a) arc also responsible for the larger shock radii 
in the 2D models G15 and M15 compared to G15-1D and M15-1D 
at early times, when multi-dimensional effects in the gain region 
do not yet play a significant role. (The data for M15-1D have been 
provided by L. Hiidepohl.) 



4. EXPLOSION DYNAMICS AND ENERGETICS 

The relativistic supernova simulations of the 11.2Mq 
and 15M0 stars discussed in this paper both yield suc- 
cessful explosions. In this section, we address the evo- 
lution of models Gil and G15 in a descriptive manner, 
focusing on the propagation of the shock, the activity 
of the SASI, the morphology of the explosion, and the 
ejecta composition. The heating conditions in our mod- 
els and the crucial differences responsible for the different 
outcome of the 15Mq models G15, M15, N15, and S15 
will be analyzed in the next section. 

Both the II.2M0 and the 15Mq models show the de- 
velopment of prompt post-shock convection a few mil- 
liseconds after bounce, which then triggers early SASI 
activity at low amplitude for « 50 ms. This can be seen 
in Figure 131 where we show the normalized coefficients 
of the decomposition of the angle-dependent shock posi- 
tion Vshid) into Legendre polynomials up to Z = 3. Ui is 
computed as 

TT 

2/4-1 r 

a^ = ^— / r,h(0)P,dcos0. (1) 


Hot-bubble convection then starts ~ 70 ... 90 ms after 
bounce and slowly pushes the shock further out than in 
ID (see Figure [SI). It again provides the seed for fast 
growth of further SASI activity, but it is only at the 
time when the Si/SiO interface reaches the shock that the 
SASI really starts to become vigorous (Figure S]). Once 
these multi-dimensional effects dominate the dynamics, 
the 11. 2 AIq and the 15Af0 progenitor evolve rather dif- 
ferently. 

4.1. 11.2A/0 Model 
4.1.1. Shock Propagation and Explosion Geometry 

In the case of the II.2M0 model, the Si/SiO-interface 
reaches the shock w 100 ms after bounce, and the as- 
sociated drop in the mass accretion rate leads to strong 
shock expansion, thus essentially enabling the approach 
to an the explosion as illustrated by Figure [5j Dur- 
ing this phase, the post-shock flow becomes dominated 
by the low-^-modes of the SASI (Figure [1]). As shown 
in Figure [3 the buoyant convective plumes that con- 
tain neutrino-heated gas merge into 2-3 large bubbles. 
Around 280 ms after bounce, the tenuous polar plumes 
disappear almost completely for a short while. This 
change of the flow geometry even results in a retraction of 
the average shock radius for 70 ms (left panel of FigurelHl 
red lines). From 350 ms onward, however, the oscillations 
of the shock become less violent as it is pushed steadily 
outward and re-sphericized somewhat by sweeping up 
mass from the spherical progenitor layers (Figures [7] and 
ini). This is also reflected by the normalized Legendre co- 
efficients ai/ao of the shock position, which decrease to 
a level of about ^ 0.1 after 450 ms post-bounce. From 
around 650 ms onward, we see positive post-shock ve- 
locities along the entire shock front. Towards the end 
of the simulations, we observe two high-entropy bub- 
bles in the northern and southern hemisphere, a rather 
broad downflow around the equatorial plane, and a fur- 
ther downflow near the south polar axis. At this stage, 
the deformation of the shock, expressed by the ratio of 
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time after bounce [s] time after bounce [s] 

Figure 4. First, second, and third coelEcient ai, a2, and as, of the decomposition of the shock radius into Legendre polynomials, 
normalized to the coefficient ao (i.e. the average shock radius) for model Gf f (left panel) and model G15 (right panel). 




0.1 0.2 0.3 0.4 0.0 0.1 0,2 0.3 0.4 0.5 0,6 

time ofter bounce [s] time after bounce [s] s[k|,/nucleon] 

Figure 5. Entropy along the north and south polar axis as a function of time for models Gil (left) and G15 (right). The trajectory of 
the Si/SiO-interface - initially located at 1200 km and « 1800 km, respectively - is visible as a discontinuity in brightness. We point 
out that model Gil (left) still shows strong SASI activity with repeated phases of shock stagnation or retraction during its long approach 
to the explosion. It should also be noted that the entropy of the buoyant bubbles pushing the shock further out in the polar directions 
is quite moderate in this case. Model G15 (right) shows a strong dipolar asymmetry during the explosion phase: While neutrino-heated 
high-entropy plumes have pushed the shock outward with a propagation velocity of > 10000 km s~^ along the northern polar direction, the 
shock radius in the opposite direction grows much more slowly and non-monotonically. Note that only parts of the simulations are covered 
by the plots because many features of the earlier phases would not be recognizable on a different scale. 




time after bounce [s] time after bounce [s] 



Figure 6. Left panel: Maximum, minimum and average shock radius for models Gil and G15. In both cases, the shock expands and 
SASI activity increases considerably when the Si/SiO-interface reaches the shock (at ^ 100 ms and fa 150 ms, respectively). It recedes 
again in the case of model G15, and the explosion is launched only several hundreds of ms later. Right panel: The "explosion energy" for 
both models (see text for exact definition). 
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Figure 7. Snapshots of the evolution of model Gil, depicting the radial velocity Vr (left half of panels) and the entropy per baryon 
(right half of panels) 115 ms, 203 ms, 290 ms, 490 ms, 658 ms, and 920 ms after bounce (from top left to bottom right). 
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Figure 8. Electron fraction (left half of panel) and density (right 
half of panel) at a time of 920 ms after bounce in model Gil. 
The white curves denote the boundary of the region where the 
local binding energy (Equation [2} is positive, i.e. they enclose the 
material that is preliminarily classified as ejecta and contains the 
high-ye {y<^ ~ 0.52) bubbles of neutrino-heated matter. 

the maximum and minimum shock radius is still appre- 
ciable: 7'sh,inax/''sh,inin ~ 4400 km/3300 km « 1.3. Even 
during the later stages of the simulation, the geometry 
of the hot plumes does not freeze out, and new down- 
flows may still develop (cp. the snapshots at 490 ms and 
658 ms). Accretion onto the proto- neutron star there- 
fore continues until late times, leading to a growth of the 
baryonic mass of the proton-neutron star from 1.275Af0 
to 1.36Af0 between 200 ms and 920 ms. 

Figure |4] (left panel) shows that both the dipole {i = 
1) and quadrupole {i = 2) mode are present in similar 
strength in model Gil, confirming that the deformation 
of the shock is largest between 200 ms and 450 ms after 
bounce and then decreases to a level of ~ 0.1 for the 
quadrupole during the subsequent evolution. 

4.1.2. Explosion Energy 

We can compute a diagnostic "explosion energy" by 
integrating over the material with positive binding en- 
ergy Cbind at a certain time. Since this energy does not 
account for subsequent nuclear recombination, burning, 
and the gravitational binding energy of the outer layers 
of the star, this quantity does not provide a direct mea- 
sure for the final supernova explosion energy. In the GR 
case, we define ebind in terms of the lapse function a, the 
rest-mass density p, the specific internal energy e, the 
pressure P and the Lorentz factor W as follows, 

ebind = a (p(c2 + e + P/p) - P) - pWc^ . (2) 

In order to mainta i n consistency with previous studies 
(iBuras et al.ll2006al: IMarek fc Jankall2009t iBnIenn et al.l 
I2009D . we exclude rest-mass contributions to the specific 
internal energy e. It can easily be verified that Equa- 
tion ([2]) correctly reduces to 

ebind P(e + P^'VS + $) (3) 

in the Newtonian limit (where $ is the gravitational po- 
tentialjB The "explosion energy" is then computed by 

^ Precisely speaking, we have a— ^l + <I>/c^ and W ^ 1 + /2 
in the Newtonian limit, and obtain the Newtonian expression as 
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integrating over the region where ebind is positive, 

£^oxpi = j ebind dF. (4) 

ebind>0 

Here, dV^ is the three-volume element for the curved 
space-time metric (and not the fiat-space volume ele- 
ment). 

The time evolution of E'expi is plotted in the right panel 
of Figure IH which shows that material behind the shock 
first becomes nominally unbound 200 ms after bounce for 
model Gil. This corresponds to the time when the shock 
first expands beyond ^ 400 km, allowing the tempera- 
ture to drop sufficiently for nucleon recombination to a- 
particles to set in. The explosion energy slowly increases 
rather unsteadily at an average rate of 6 x 10^^ erg s^^, 
and then seems to level off around 3.5 x 10"*^ erg around 
simulation after 600 ms post-bounce with some residual 
fluctuations. By the end of simulation, the total mass 
of the material with positive binding energy amounts to 
O.O45M0. 

Despite this seeming "saturation" of ii'cxpi, no defini- 
tive statement about the final explosion energy can be 
made as yet at this stage, although we can tentatively 
estimate corrections due to the energy input from nu- 
clear burning in the shock, £^burn, the recombination 
of nucleons into a-particles (and possibly further into 
heavy nuclei), and the deduction of the binding energy 
^'proshock of the uushockcd matter. All of these cor- 
rection terms are of a similar magnitude as the cur- 
rent diagnostic value i?expi: -E'preshock yields a contribu- 
tion of roughly —7.5 x 10^^ erg (i.e. exceeding i?expi), 
but this value is beset with uncertainties as part of the 
shocked material is still accreted onto the proto-neutron 
star (which will reduce the energy required to unbind 
the outer shells). Recombination of nucleons and a- 
particles in the electa would provide an additional energy 
of E-ccc « 2 x 10* erg. Burning in the shock will not yield 
any significant contribution with the current shock veloc- 
ities as the typical post-shock temperatures are already 
too low (< 3 X 10^ K) to allow for explosive Si- and O- 
burning. The uncertainties in these numbers illustrate 
that in order to obtain a reasonably accurate prediction 
for the explosion energy of model Gil, the simulation 
probably needs to be extended until the shock reaches 
the C-0 shell at > 2 x 10* km (at which point the bind- 
ing energy of the remaining pre-shock matter would be 
negligible) . The final explosion energy depends critically 
on the fraction of shocked material from this shell that is 
accreted onto the proto-neutron star and thus need not 
be unbound at all. 

The slow growth and the stagnation of the explosion 
energy visible in Figure [6] at late times is a consequence 
of relatively inefficient neutrino heating in model Gil 
associated with an unfavorable 2D fiow geometry in this 
case. In the late phases, the rate of energy input by neu- 
trinos in the gain region, Q^, is only '-^ 5 x 10^° erg s~^. 
With a typical binding energy of a mass element at the 
gain radius of ~ 30 MeV/baryon, this implies that only 
^ O.OIMq s~^ of additional material can be unbound by 
neutrino heating. The actual mass fiux from the heating 
region into the ejecta is somewhat higher because re- 

an approximation to 0{e/ , P/ pc^ ,v'^ / / c^). 
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combination of nucleons into a-particles also contributes 
part of the energy for unbinding the newly ejected mate- 
rial. Since the ejecta from the gain region are channeled 
through relatively narrow outflows into high-entropy 
bubbles (cp. Figure [7]) at high velocities, neutrino heat- 
ing is rather inefficient due to the short exposure time. 
This also implies that only a small excess energy - i.e. 
much smaller than the maximum energy available from 
recombination of 8.8 MeV/baryon - remains for increas- 
ing the total (internal-|-kinetic~|-gravitational) energy of 
the ejecta. Newly ejected material thus adds only a few 
10**^ erg s~^ to the explosion energy. However, the shock 
also sweeps up material at a rate of about COSM© s~^, 
which implies a negative energy flux into the "ejecta re- 
gion" of 5 X 10*^ erg s~^. This may balance or even 
exceed the energy carried by fresh ejecta from the gain re- 
gion, thus accounting for the unsteady evolution of -Ecxpi- 
With energy being fed into the ejecta at such a low rate, 
a considerable fraction of the material swept up by the 
shock will be channeled into the downflows and accreted 
onto the proto- neutron star (Figure ^ . We already ob- 
serve that the accretion rate onto the proto-neutron star 
starts to increase again towards the end of the simulation, 
which leads to a late-time rise of the electron neutrino 
and antineutrino luminosities (which will be discussed in 
a subsequent paper on the neutrino signal). All these in- 
dications suggest that model Gil remains a low-energy 
case and is likely to represent a fallback supernova. 

4.1.3. Ejecta Composition 

Although the final ejecta composition for our mod- 
els can only be determined by detailed nucleosynthe- 
sis calculations once the amount of fallback is known, 
a few remarks can already be made about the nucle- 
osynthesis conditions on the basis of the entropy (Fig- 
ure [7]) and the electron fraction (Figure [8]) of the ma- 
terial that is presumably ejected. Low-entropy matter 
(s < 10 kf,/nucleon) with ~ 0.5 that has under- 
gone little or no neutrino heating, but has either been 
swept up directly by the shock or has been pushed out by 
neutrino-heated bubbles before falling inward to smaller 
radii, contributes most of the mass in the "ejecta region" . 
Depending on the maximum temperature reached be- 
fore expansion sets in, this material has partially been 
processed by nuclear burning to various degrees: Ac- 
cording to the simple "flashing" treatment in the Ver- 
tex code (Ranipp & Janka 2002) , intermediate- mass el- 
ements dominate the composition. 

Hot, neutrino-processed material with entropies of up 
to 35 kf,/nucleon makes up only for a small fraction of 
the ejecta (- 0.005Mq or ~ 11%). This part of the 
ejecta is proton-rich with an electron fraction Y,, rang- 
ing from sa 0.51 up to w 0.58. Different from the 
case of elec tron-capture super novae (Wanaio et al.ll201ll ) 
and unlike iPruet et"al] ()2005[) , who considered an artifi- 
cially triggered 2D explos ion of the 15M0 progenitor of 
I Wooslev fc Weaved (flOOSh . we do not observe any slightly 
neutron-rich pockets in the ejecta, which is a consequence 
of the different (slower) outflow dynamics (a detailed 
analysis will be provided in a forthcoming paper). We 
therefore expect the nucleosynthesis yields to conform 
with the established results for proton-rich outflow con- 
ditions, i.e. the ejecta composition will be dominated by 
^^Ni and helium with an admixture of a few rare iso- 



topes (^fSc,_^flj_,a£d ^^Zn)_w^ fac- 
tors ()Pruet et al.ll2005t iFrohlich et al.ll2006bl) . Depend- 
ing on the neutrino luminosities there m ay also be poten- 
tial for t/jp-process nucleosynthesis fjFrohlich et al.l[2"006al : 
IPruet et al.l[2006h . 

4.2. 15Mo Model 
4.2.1. Shock Propagation and Explosion Morphology 

Different from the 11.2M0 case, the development of the 
explosion in model G15 is not immediately connected to 
the Si/SiO interface reaching the shock although the de- 
crease in the accretion rate results in a transient increase 
of the average shock to 220 km at 200 ms after bounce 
(see Figure |3] and Figure [5]) and also in increased SASI 
activity (Figure HI right panel). However, 230 ms after 
bounce the shock again starts to recede slowly with the 
average shock radius reaching a minimum value of about 
100 km at 380 ms. During this period, the SASI remains 
active with strong dipole and quadrupole components 
(the maximum amplitudes being ai/ao « 02/00 ~ 0.3, 
FigureSl right panel). Around 400 ms, the average shock 
radius begins to move outward rather steadily (Figure[3l), 
and at about 430 ms, some material becomes nominally 
unbound (Figure [6]). Model G15 develops a strongly 
asymmetric explosion (Figures I5|6|9p : By the end of the 
simulation, the shock has reached 3800 km in the north- 
ern hemisphere, while the minimum shock radius over the 
only remaining strong downflow in the southern hemi- 
sphere is only 850 km (Figure |6|); i.e the ratio r-max/'^min 
of the maximum and minimum shock radius is as large 
as 4.5:1. Snapshots of the developing asymmetric explo- 
sion with even more extreme shock deformation during 
earlier phases of the explosion are shown in Figure [HI 

4.2.2. Explosion Energy 

From 540 ms onward, we observe a much more rapid 
growth of the explosion energy than for model Gil. By 
770 ms, it has reached a value of i^oxpi = 1.3 x 10^" and 
still continues to increase at a rate of w 7 x 10^° erg s~^. 
This rate of increase is roughly consistent with the as- 
sumption that neutrino heating at a rate of (4 ... 5) x 
10^^ erg s"^ (Figure 1131) allows an outflow rate from the 
gain region of ~ O.O6M0 (assuming a binding energy 
of ^ 40 MeV/baryon at the gain radius), and that recom- 
bination of nucleons into a-particles leads to an excess 
energy of ~ 8 MeV/baryon that actually contributes to 
the explosion energy. As for model Gil, more material 
is actually swept up by the shock (~ O.15M0 s~^) than 
is supplied from the heating region, but in contrast to 
model Gil, the corresponding negative total energy flux 
into the ejecta region is much smaller (~ 10^" erg s~^) 
than the positive contribution from hot, neutrino- heated 
ejecta (~ 10^"'^ erg s~^). -Ecxpi therefore stably increases 
for model G15. The presence of a strong and stable 
downflow in the southern hemisphere is probably help- 
ful for this behavior: Compared to model Gil, the ac- 
cretion rate onto the proto-neutron star is considerably 
higher (M ~ 0.1 M0 s~^) during the late phases, which 
results in a sizable accretion contribution to the neutrino 
luminosity and therefore allows for persistently strong 
neutrino heating in the gain region. Due to continuous 
accretion, the baryonic mass of the PNS has grown to 
I.58M0 by the end of the simulation. 
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Figure 9. Snapshots of the evolution of model G15, depicting the radial velocity v,- (left half of panels) and the entropy per baryon 
(right half of panels) 135 ms, 207 ms, 332 ms, 534 ms, 644 ms, and 775 ms after bounce (from top left to bottom right). 
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Figure 10. Electron fraction (left half of panel) and density (right 
half of panel) at a time of 775 ms after bounce for model G15. 
The white curves denote the boundary of the region where the 
local binding energy (Equation [2]| is positive, i.e. they enclose 
the material that is preliminarily unbound and classified as ejecta 
(mostly located in the northern hemisphere) and contains the high- 
Ye {Ye ~ 0.52) bubbles of neutrino-heated matter. 

With respect to the final explosion energy, the same 
reservations apply as for model Gil, i.e. the binding en- 
ergy of the pre-shock material (roughly 2.6 x 10^*^ erg) is 
still larger than the current value of £^oxpi- However, 
different from model Gil we expect a significant en- 
ergy input from burning in the shock on the order of 
(1 ... 2) X 10^° erg or more, i.e. model G15 has almost 
reached the stage where the envelope could be unbound 
completely even without further neutrino heating. 

4.2.3. Ejecta Composition 

The composition of the ejecta in model G15 exhibits 
some marked differences to model Gil, although we can 
identify the same two components: low-entropy isospin 
symmetric matter swept up by the shock, and proton- 
rich, neutrino- heated matter (Figures I9ll0p . Different 
from model Gil, much of the shocked material accu- 
mulated by the large expanding bubble in the northern 
hemisphere has been burnt to ^^Ni, and even at 775 ms 
post-bounce, the post-shock temperature is still high 
enough to allow at least for explosive 0-burning. More- 
over, hot, neutrino-heated material accounts for a larger 
fraction of the ejecta (O.OIMq, or roughly one third of 
the matter with positive total energy). We again find 
the neutrino-heated ejecta to be exclusively proton-rich, 
but both the typical and the maximum electron fraction 
in the neutrino-heated ejecta are even higher than for 
model Gll with Y,, ranging up to 0.6. 

5. ANALYSIS OF HEATING CONDITIONS AND 
COMPARISON TO NON-EXPLODING MODELS 

It is interesting to compare the evolution of model G15 
to the pseudo-Newtonian run M15, the purely Newtonian 
model N15, and to simulation S15 with simplified neu- 
trino reaction rates. The average shock radius is rela- 
tively similar in all cases until « 400 ms after bounce 
(Figure except for model N15, which maintains a 
larger average shock radius (by 30 ... 50 km) . Among 
G15, M15, and S15, the shock in model G15 reacts most 
strongly to the drop in M associated with the Si/SiO 
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interface and stays a little further out until 380 ms, but 
the differences between the three models remain rather 
modest; right before the onset of the explosion in model 
G15, the shock position is virtually identical. Yet, only 
a few tens of milliseconds later, an explosion develops 
in model G15, whereas S15, M15 and even N15 with its 
relatively large shock radius (« 140 km) show no sign 
of shock revival. In this section, we discuss the reason 
for this (unexpectedly) different behavior of the models 
lacking either the GR hydrodynamics treatment or the 
full set of neutrino opacities. 

5.1. Analysis Framework 

Shock revival can be understood in terms of 
a competition of neutrino heating and the down- 
ward advection of the gas in th e gain layer (e.g. 
iBurrows fc GosiwHTOOl lJankal[200l . A detailed anal- 
ysis shows that there is a critical neutrino luminos- 
ity for each value of the mass accretion rate M of 
the shock, above which no steady-state accretion so- 
lution can exist ([B urrows & Goshv 1993; Janka 200fl 
'Murphv fc Bur rowi 120081: iPeicha fc Thompson! I2012t 
iFernandez 20li[ ). Roughlv speaking, conditions become 
favorable for an explosion if the time the accreted ma- 
terial spends in the gain region (measured by the "ad- 
vection", "residency", or "dwelling" time-scale Tadv) be- 
comes longer than the time required to lift its binding 
energy to positive values by neutrino heating (the heat- 
ing time-scale Thoat)- If a ratio Tadv/Thcat ?J 1 is main- 
tained for a sufficiently long time (typically a few tens 
of milliseconds), the shock can, in all probability, ex- 
pand sufBciently to create a positive feedback loop by 
further increasing Tadv, thus establishing a runaway sit- 
uation with sustained shock expansion that eventually 
leads to an explosion. This concept provides an ade- 
quate basis for understanding the different outcome of 
the exploding models Gll and G15 in con trast to the 
other 15AfQ runs fcp. IMarek fc Jankal [20091 ). 

We compute Tadv as the ratio of the binding energy 
I-Egaini of the material in the gain region and the volume- 
integrated energy deposition rate in that region (in agree- 
ment with IMarek fc Jankal 120091 and with the best defi- 
nition identified bv lFernandedl2011l ). 



J heat 



(5) 



Here, i?gain and Qhcat are volume-integrals over the bind- 
ing energy density Cbind as given by Equation ^ and 
the local neutrino heating rate per unit volume be- 
tween the gain radius Tgain and the (average) shock ra- 
dius Tshock, 



GbinddT^, 



ai„<I'<r3i,ock 



; heat 



(6) 



(7) 



'"gain<''<''sho 



There are alternative definitions for the advection time- 



scale Tad,, (iBuras et al.l 


2006al IMarek & Jankal 120091: 


iMurphv fc Burrowsl 12008 


; Peicha fc ThompsonI l2012f). 



each of which can be supported by very plausible argu- 
ments. For the sake of convenience, we choose to express 



Relativistic Explosion Models of Core-Collapse Supernovae 



13 



it in terms of two readily available quantities, the accre- 
tion rate M of gas through the shock and the baryonic 
rest mass contained in the gain layer Mgain, 

= Afgain/Af, (8) 

where A/gain is given by 

A//gain= j pWdV. (9) 

'•gain<i'<rshock 

This definition measures the time matter needs to flow 
through the gain region if steady-state conditions hold, 
and is used here because it allows for a very straightfor- 
ward evaluation of Tadv/'^'hcat not only in the ID case, but 
also in the multi-dimensional case. The time-scale ratio 
■'"adv/'^'hoat was recently also shown by iFern andcz (20lT|) 
to provide a useful instrument for distinguishing models 
that are going to explode from "pessimistic" ones. 

5.2. 11.2A^Q Progenitor 

The evaluation of the time-scale ratio Tadv/'''hcat con- 
firms that the transition to shock expansion of model Gil 
is a relatively clear-cut case, and is indeed triggered by 
the Si/SiO-interface. As shown by Figure[TTl the expan- 
sion of the shock due to the drop in the mass accretion 
rate results in a considerable increase of the advection 
time-scale Tadv (by a factor of 2 ... 3) within a few tens of 
milliseconds, which is sufficient to bring the critical time- 
scale ratio Tadv/Thcat above unity in conjunction with 
the support from hot-bubble convection. At later times, 
Tadv/^heat uevcr drops below unity, not even during the 
transient phase of shock retraction around 300 ms, in- 
dicating heating conditions that favor robust shock ex- 
pansion. These findings are in v ery good agreement with 
the pseudo-New t onian model of lBuras et al.l (|2006al ) and 
IMarek fc Janfal (1200^ . even though the explosion mor- 
phology is different (spherical instead of dipolar), which 
is probably the result of stochastic variations of the SASI 
activity. 

5.3. 15Af0 Progenitor ~ Influence of General Relativity 
and Neutrino Interaction Rates 

The case of the 15Mq model is much more interesting 
because the different outcome of model G15 as opposed 
to M15, N15, and S15 suggests that this progenitor is 
a marginal case and may reveal the dependence of the 
explosion conditions on the input physics used in these 
three runs, i.e. the influence of the GR treatment and 
the neutrino reaction rates. Indeed, a comparison of the 
runaway criterion Tadv/'''hoat and the critical time-scales 
(Figure fTTj) shows systematic differences: For model G15, 
we consistently find higher values of Tadv/^heat than in 
the other 15Mq runs, and this model is also charac- 
terized by a larger heating efficiency Q\ica.tl {Lu^ + LuS) 
(Figure fT2|) . i.e. a larger fraction of energy radiated in 
Ve and Ve is re-absorbed in the gain layero Moreover, 
the volume-integrated heating rate Qhcat is also highest 
among the 15Af0 runs (Figure [T3l) . 

Compared to model M15, Tadv/Theat is typically higher 
by a factor of two, resulting from the combination of 

^ Note that in contrast to IMarek fc Jankal II2009I ). we use the 
luminosities of v,^ and Pe at the gain radius (see Figure [l5]l instead 
of the (lower) values for an observer at infinity. 
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Figure 11. Runaway criterion T^dv /'''heat i heating time-scale 
Thcati Enid advection time-scale r^dv for models G15 (thick black 
solid line), M15 (thick, red, dashed), S15 (thick, blue, dash- 
dotted), N15 (thick, black, dotted), and Gil (thin solid line, black) 
smoothed over 10 ms. Note that we do not evaluate the time-scales 
during the later phase of the explosion when the their definitions 
are no longer meaningful. Naturally, the time-scales cannot be 
computed heiore the formation of the gain layer several tens of 
milliseconds after bounce either. 
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Figure 12. Heating efficiency for models G15 (black solid line), 
M15 (red, dashed), S15 (blue, dash-dotted), and N15 (black, dot- 
ted). The heating efficiency is computed as the ratio '3hcat/(^i^e + 
Li?^ ) of the volume-integrated neutrino heating rate Qheat S'^d the 
sum of the Ug and luminosities. 

both a longer advection time-scale and a shorter heating 
time-scale. The purely Newtonian model N15 exhibits a 
very similar evolution of the runaway criterion as model 
M15, but this is due to a cancellation of huge differences 
(amounting to a factor of 2 ... 3) in both the advection 
and the heating time-scale. The relativistic calculation 
shows, however, that even though the effects of a smaller 
shock radius (Figure [3]) and hence a shorter advection 
time-scale on the one hand and stronger heating on the 
other hand partially compensate each other in GR, such 
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Figure 13. Tlie energy deposition rate in the gain region (top left panel, smoothed over 10 ms), the total (time-integrated) energy 
deposited in the gain region (bottom left), the total mass in the gain region (top right, smoothed over 10 ms), the gain radius and the 
maximum shock radius (fjottom right) for models model G15 (black solid line), M15 (red, dashed), S15 (blue, dash-dotted), and N15 (black, 
dotted). 
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Figure 14. Circumferential radius (decreasing with time) and sur- 
face temperature (increasing) of the proto-neutron star for models 
G15, M15 S15, and N15. The surface is defined by a fiducial density 



a cancellation is not to be taken for granted, and the 
residual effect, e.g. on the time-scale ratio Tadv/Thcat niay 
still be on the order of several tens of percents. 

For model S15, the time-scale ratio lies about halfway 
in between G15 and N15, suggesting a non-negligible ef- 
fe ct of the neutri n o mic rophysics. Considering the claims 
of IBruenn et al.l ([2009) to that effect, this is certainly 
noteworthy, but we shall first turn our attention to the 
more pronounced differences between the GR model G15 
and the runs M15 and N15 with a different treatment of 
gravity. 

5.3.1. Newtonian Approximation vs. General Relativity 

The purely Newtonian case (N15) stands apart most 
clearly from the others with large values of Tadv and 
Theat, and the reason for this has essentially been given 



bv IBruenn et al.l (|2001[ ): Due to the shallower potential, 
the proto-neutron star is considerably more extended in 
model N15 compared to the GR case (Figure [HI) , and 
this also shifts the gain radius and shock radius further 
out (Figures |3] and [T3l) , thus increasing the mass Mgain 
in the gain layer. On the other hand, the neutron star 
surface temperature is also considerably lower compared 
to GR (Figure [T4l) , resulting in a significant reduction 
of the neutrino luminosities and mean energies in the 
gain region (Figure [T5)) and hence weaker heating in the 
gain region (Figure [TS]). Apparently this also leads to 
less vigorous convection in the purely Newtonian case as 
the kinetic energy contained in non-radial mass motions, 
Ekin,e, is typically lower than in model G15 (Figure [T6)) . 
The relative weakness of convection in model N15 is in- 
dicated even more clearly by the mass-specific kinetic 
energy £'kin,e/Afgain, which is lower by 40% — 50% than 
in model G15. Stronger convection in GR partly com- 
pensates for the reduction of Mgain and Tadv due to the 
smaller neutron star radius and therefore helps to turn 
the scales in favor of a larger value of the runaway crite- 
rion Tadv/Thcat " the large reduction of Thoat emerges as 
the dominant effect (Figure [Tl]) . 

In comparing model N15 to G15 and M15, we should 
also bear in mind that N15 was computed with a higher 
angular resolution of 128 zones, which might be beneficial 
for the heating conditions bec ause it was seen to foster 
explosions in 2D simulations of i Hanke et al.l (|201lD . The 
increase of Tadv/Theat in GR compared to the Newtonian 
approximation may therefore even be underestimated by 
our analysis. 
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Figure 15. Neutrino mean energies (upper panels) and luminosities (lower panels) at the gain radius for the relativistic model G15 
(black solid lines), the pseudo-Newtonian model M15 (red dashed lines), the relativistic model S15 with simplified neutrino rates (blue 
dash-dotted lines), and the purely Newtonian run N15 (black, dotted). Electron neutrinos, electron antincutrinos, and /x/t neutrinos are 
shown in the left, middle, and right panel of each row, respectively. The mean energy is defined as the ratio of the angle-averaged neutrino 
energy and number density, and the luminosity is computed as the integral of the fiux over the sphere corresponding to the gain radius. 
The luminosity given here is for a single species of i^^/ry for all four of them. 

' ' ' ' ' ' " ' Interestingly, the tendency towards slightly more ener- 

getic z/o's and D^s in GR is already present in ID (see Pa- 
per I) . This is presumably the result of a slightly different 
density stratification in GR that cannot be reproduced 
exactly by the modified Newtonian potential and the ap- 
proximate GR transport treatment (e.g. due to the iden- 
tification of coordin ate radius and pro per radius) used 
for the M15 run (cp. lMarek et al.ll2006[ ). The circumfer- 
ential radius of the PNS (defined as the radius where the 
density drops to 10^^ g cm~'^) is indeed larger by 2%-4% 
in GR, and its surface is somehwat hotter at late times 
(Figure [T4|) . In contrast to the higher i/^ and De luminosi- 
ties in model G15, the luminosity of fi and t neutrinos is 
smaller in the GR case. This is the result of general rel- 
ativistic transport effects and a different stratification at 
the edge of the PNS convection zone, where the "num- 
ber sphere" for and i>t (i-e. the sphere up to which 
production processes and abs orption processe s remain in 
equilibrium; seelRaffeh 20011: iKeil et~aIll2003D is located, 
as we illustrate for a representative snapshot 250 ms after 
bounce (Figure fTT]). The upper panel of Figure [TTl shows 
the energy fiux density of i^^/r in the comoving frame 
for models G15 and M15 at that time, with H^, rising to 
somewhat higher values before it decays like outside 
the "number sphere" To demonstrate that the smaller 
neutrino flux in the case of model G15 is consistent with 
the temperature and density stratification depicted in the 
lower panel of Figure [17] we also show a rough estimate 
for the diffusive fiux for comparison in the upper panel. 

* We remark that the luminosity emerging from the "number 
sphere" is smaller than the black body luminosity by one to two 
or ders of magn itude due to a very small flux factor (see discussion 
in IJankal 1 19951 1. Stefan's law is therefore ill-suited for estimating 
the luminosity of i/^ and Ut- 
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Figure 16. The kinetic energy -Ekin.e associated with lateral mo- 
tions in the gain region for models G15 (black solid line), M15 
(red, dashed), S15 (blue, dash-dotted), and N15 (black, dotted); 
the curves are smoothed over 10 ms. i?kin,e can be taken as a mea- 
sure for the combined activity of convection and the SASI, which 
both contribute to this quantity. 



5.3.2. 



Effective Potential Approximation vs. 

Relativity 



General 



The comparison between models G15 (GR) and M15 
(effective potential) is somewhat more subtle, but the 
different heating conditions can still be traced back - at 
least partly - to the neutrino emission from the proto- 
neutron star. Again, the neutrino luminosities and mean 
energies at the gain radius (Figure [15]) turn out to be the 
crucial factor. In the early phase, the GR run exhibits a 
noticeable enhancement in the electron antineutrino lu- 
minosity (by ~ 15%) and to a lesser extent in the electron 
neutrino luminosity. In addition, the mean energies of Vc 
and Pe tend to increase more strongly at late times, with 
the difference reaching almost 1 MeV for the antincutri- 
nos. 
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Figure 17. Origin of the different and luminosities in mod- 
els G15 and M15. The upper panel shows the energy flux density 
of and Vt at 250 ms after bounce (black lines) for G15 (solid) 
and M15 (dashed) and illustrates that a smaller maximum flux is 
reached at a somewhat smaller "freeze-out" or "number sphere" ra- 
dius in model G15 (corresponding to a saturation of the luminosity 
at a lower value). Dotted lines are used to denote the rough loca- 
tion of the freeze-out radius in both models. The lower flux at radii 
around 24 km is consistent with an estima te o f the flux H^^^j^ ex- 
pected for equilibrium diffusion (Equation IIOI I as indicated by the 
red curves: While the diffusive fluxes agree well within a large part 
of the convection zone, the shift of the "dip" in the predicted //^ dift 
in G15 suggests a lower freeze-out radius and luminosity. Note that 
the diffusive flux has been rescaled to avoid overlap with the black 
curve. Since the density- and temperature-dependence of the dif- 
fusion coefficient cannot be captured exactly by simpl e po wer laws, 
a perfect fit is not to be expected anyway; Equation mOI I is rather 
used to detect relative differences between the two models G15 
and M15. An gle- averaged profiles of the quantities used to evalu- 
ate Equation mOI I. i.e. the density p, the temperature T, and V are 
shown in the lower panel along with the angle-averaged entropy 
s. The different fluxes in models G15 and M15 stem from several 
(partially competing) factors, namely a shallower temperature gra- 
dient in model G15 from 22 km outward, which is associated with 
an earlier rise of s outside the proto-neutron star convection zone, 
lower densities (which increase the diffusive flux; the relative dif- 
ference is still ~ 10% at 25 km), and the reduction of the flux due 
the F-factor (which is not present in M15). The lower freeze-out 
radius in model G15 is associated with a different location of the 
edge of the (roughly adiabatically stratified) PNS convection zone, 
which is clearly visible in the entropy profile. 

The diffusive fltp c i?,ydiff can be computed as (see, e.g., 
iPons et al.|[T99l Equation 130 



,diff 



kD- 



y3 



(10) 



using Schwarzschild radial coordinates (with r denoting 
the circumferential radius) . Here a is the lapse function, 
D an appropriate, frequency-averaged diffusion coeffi- 
cient, and some factors of order unity have been absorbed 

^ In Vertex-CoCoNuT, we do not distinguish between heavy 
fiavor neutrinos and antineutrinos, and the chemical potential 
P^v^l^ therefore vanishes. Consequently, only derivatives of the 
temperature appear in the diffusion equation. 



into the constant fc for convenience. A metric factor F 
appears because the diffusive flux depends on the deriva- 
tive Td/ dr with respect to proper (physical) radius^ In 
model M15, no distinction is made between proper radius 
and circumferential radius (F = 1). As scattering on nu- 
cleons is the dominant source of opacity, the diffusion 
coefficient D is roughly proportional to p^^T"^ in the 
density regime that we are considering here. Figure [17] 
demonstrates that Eq. (|10p predicts different fluxes in the 
region around the "number sphere" and also points to a 
slightly smaller "number sphere" radius in G15, which is 
associated with a different location of the outer boundary 
of the PNS convection zone (Figure [T7|) . The lower flux, 
mainly due to the factor F in the GR case and the smaller 
"number sphere" radius account for the ~ 20% lower 
and h'r luminosity in model G15. On the other hand, 
there is no such reduction of the neutrino mean energy of 
and Vr in model G15, because neutrino-electron scat- 
tering and non-isoenergetic neutrino-nucleon scattering 
still allow for energy exchange and thermal equilibration 
with the medium at larger radii, where the temperature 
is higher than in model M15. 

While GR transport and stratification effects thus seem 
to account for the different neutrino emission in mod- 
els G15 and M15 and, because of the higher and 
De luminosities, suggest somewhat better heating con- 
ditions, the huge increase of the time-scale criterion in 
model G15 still needs to be connected to these find- 
ings. Merely reducing the heating time-scale by « 20% 
in model G15 compared to M15 would not be sufficient 
to ensure that the runaway condition Tadv/Thcat > 1 is 
met. At 450 ms the crucial factor distinguishing the rela- 
tivistic and the pseudo-Newtonian run is the larger value 
of Mgain (Figure [131) and the longer advection time-scale 
T^adv = Mgi^in/M (Figure [TT|) . Part of this increase in 
Tadv is a direct consequence of the higher neutrino lu- 
minosities, as the increase in thermal pressure behind 
the shock allows for a larger mass in the gain region. 
Since the position of the gain radius is almost identical 
in model G15 and M15 (Figure [T5)) . this can be illustrated 
qualitatively on the basis of a very rough approximation: 
Requiring balance between heating and cooling at the 
gain radius with heating and cooling rates per nucleon 
roug hly proportional to L(Ef,) /r"^ and T^, respectively 
fe.g. IBethe fc WilsonI [1981 IJankal[200l . and assuming 
that the gas at the gain radius is radiation-dominated, 
we find that the pressure Pgain and the temperature Tgain 
at the gain radius scale as 



r)3/2 ry-iG 
-'gain ^ -'gain ^ 



L. {El) 



(11) 



The stratification in the gain region is roughly adia- 
batic with the pressure and density following power laws 
(P cx r~^, p oc r~^), and in spherical symmetry we can 
therefore approximatively determine the radius of the 
stagnant accretion shock with the help of the jump con- 
ditions at the shock. For the post-shock pressure Ppost, 

^ In Vertex-CoCoNuT we work with a different gauge choice 
for the radial coordinate (isotropic coordinates riso), but the con- 
version to the familiar Schwarzschild form of the metric is straight- 
forward, r can be obtained from the conformal factor </> as 
r = l-h2ri,oain<^/ari,o. 
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we have 



^post — ( 1 ^ ^ I 



(12) 



where Pram is the rara pressure ahead of the shock, and 
P is the ratio of the post- and the pre-shock density. 

— Ppro^^pro '^^^ be computed assuming that the pre- 
shock velocity Wpro is a certain fraction of the free-fall 



velocity, i.e. 



OC 1/, 



and therefore scales as 



-^lam - . 2 ^ 5/2 ' 



(13) 



On the other hand, the post-shock pressure is related to 
the pressure at the gain radius roughly by 

Ppost^sh ~ ^gain^'gain- (14) 

Assuming rgain to be fixed (motivated by Figure [T5)) . the 
shock radius therefore varies with L^{El) as 

rsh OC (L, {El)f^ . (15) 

Using Equation (|15p . we can compute the mass Mgain in 
the gain region for a density stratification with p cx r~^, 



Mgain OC 



I3M /r,h 



^ T y \ ^gain 



The logarithmic derivative, 

SlnAfgain 2 , 



ain(L,(S2)) 3 91n(r,h/ 



(16) 

1, (17) 



can be taken as a measure for the sensitivity of A/gain to 
changes in L^, {Ef,). Although based on a rather crude 
approximation. Equation ()17p suggests that changes in 
the luminosity and neutrino energy induce comparably 
large changes in Afgain and hence Tadv 

In addition, stronger neutrino heating also leads to 
more violent activity of non-radial mass motions in the 
form of convection and the SASI (see Figure [TO]) . This 
in turn increases the residence time-scale of matter in 
the gain region and thus further boosts the effect of the 
stronger neutrino heating in model G15. With both 
mechanisms working in combination, the mass in the 
gain region Afgain and hence the advection time-scale 
Tadv reach considerably larger values in model G15 (ex- 
cept for a short quiet period (see Figures [5l5|) with little 
SASI activity between 350 ms and 400 ms). Likewise, 
the energy deposition rate in the gain region (Figure I13p 
remains higher than in model M15 by a factor of ~ 2 
from 200 ms onward. 

While we view all these aspects as a likely explana- 
tion, it cannot be excluded that other factors are respon- 
sible for the different evolution of the shock after the 
accretion of the Si/SiO and the relatively quiet period 
around 350 ms. Different growth and saturation proper- 
ties of the SASI due to general relativistic effects could 
be invoked as an alternative explanation, although this 
seems rather far-fetched considering that the non-linear 
behavior of the SASI is still under investigation even in 
the Newtonian case (where parasitic Kelvin-Helmholtz 
or Ray leigh- Taylor instab ilities have been p roposed as a 
saturation mechanism by iGuilet et al.ll2010l) . The possi- 
bility of numerical effects (in a broad sense) ought to 



be considered more seriously: The hydro solvers Co- 
CoNuT and Prometheus used for model G15 and M15, 
respectively (Table[T]), employ different Riemann solvers, 
Prometheus relies on dimensional splitting, the radial 
grids arc not identical, and neither are the initial per- 
turbations. However, there is little grounds for relating 
these technical differences to the behavior of models G15 
and M15. Apart from the fact that the Si/SiO interface 
is preserved a little more sharply in CoCoNuT (which 
may account for a slightly stronger reaction of the shock 
to the reduced accretion rate. Figure [3l), we have no tan- 
gible evidence for such a connection. Moreover, the fact 
that the CoCoNuT model S15 - with a less pessimistic 
ratio Tadv /''"heat (Figure [Tl]) - also responds less vigor- 
ously to the drop in the accretion rate than model G15 
points to a physical origin as well. It therefore seems 
that the more rigorous GR treatment is ultimately the 
reason for the more optimistic evolution of model G15 
compared to the pseudo-Newtonian simulation M15. 

5.3.3. Influence of Neutrino Rates 

The more optimistic evolution of model G15 compared 
to the S15 run with a simplified set of neutrino interac- 
tion rates can be understood in a similar vein by consider- 
ing the neutrino emission from the PNS (Figure [T5|) . The 
most conspicuous effect in model G15 is the reduction 
of the v^/t mean energies (because of non-conservative 
neutrino- nucleon scatterings, see below), which is of lit- 
tle direct relevance for the energy deposition in the gain 
layer to which and v-r hardly contribute. However, 
electron neutrinos and antineutrinos are also affected on 
a smaller scale with some repercussions on the heating in 
the gain layer. In particular, electron antineutrinos are 
emitted with higher luminosity (by about 10%) and - to- 
wards the later phases - higher mean energies (by up to 
0.8 MeV) in model G15. For the electron neutrinos, the 
luminosity enhancement is smaller, and the spectra are 
only a little harder with the improved set of interaction 
rates. 

However, Figure [TT] shows that using the improved 
opacities results only in a relatively insignificant decrease 
of the heating time-scale in G15. This small effect is mag- 
nified by the increase of the advection time-scale due to 
a larger average s hock r adius (see Figure [3] and the dis- 
cussion in Section [5.3.2[) . as well as by increased convec- 
tive and SASI activity (see Figure [TO)) , and therefore still 
leads to an appreciable difference in the time-scale ratio 
Tadv/Tiicat between models G15 and S15. The reduced 
heating in model S15 thus at least results in a consider- 
able delay of a possible explosion in this model. 

It remains to be discussed how the improved interac- 
tion rates in model G15 contribute to an enhancement of 
Vf, and Df, emission. Upon closer analysis, energy trans- 
fer from /i and t neutrinos to nucleons by scattering re- 
actions in the PNS surface region (not in the gain re- 
gion!) emerges as the most likely ca use. As pointed out 
by several authors (Janka & Hilleb randtl 119891 : iSuzukil 
IT990t IRaffell]l200lt iKeil et al.ii2003l) . and Ur pass an 
extended scattering atmosphere (dominated by elastic 
neutrino- nucleon scattering reactions), which separates 
the "number sphere" where the production reactions (in 
our case mainly bremsstrahlung and neutrino-neutrino 
pair conversion) freeze out, and the "transport sphere" 
which marks the transition to free streaming. In this 
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Figure 18. Illustration of enhanced and Pe emission as an indirect consequence of non-isoenergetic nucleon recoil effects in the !^^/r 
sector. The two panels show the redshifted lab-frame luminosities of (left panel, black curves) and Pe (right panel, red curves) at a time 
of 400 ms after bounce as a function of the average density p(r) at a given radius r for model G15 (black solid lines, non-isoenergetic nucleon 
recoil included) and model S15 (dash-dotted lines, without energy losses of neutrinos in nucleon scatteringsq). At p(r) fs! 5 X 10^^ g cm~^ 
(indicated by the dotted vertical line), the Ve (and, respectively, the Pe) luminosities with and without nucleon recoil are roughly identical. 
For the heavy-flavor neutrinos (y^i,., Pp, Vt^ and Pr), we show halj the total luminosity of the four species and shift the resulting curves to a 
value of 3 X 10^^ erg s~^ at p(r) 5 X 10^^ g cm""^ (blue curves) in order to illustrate the different amount of energy transfer from f^^/r 
to the medium below this density. Energy transfer from i^^/^ (reflected by the decrease of L,, towards lower densities) to the background 
medium occurs roughly between 5 X 10^ g cm~^ and a few 10^^ g cm"''. With non-isoenergetic nucleon recoil, the energy input into the 
medium is larger, and there is a corresponding enhancement of the Ue and Pe luminosity in precisely the same region. If the additional 
energy gained from i'^/t scattering were completely transferred into and equally split between Ue and Pe, the expected enhancement of 
the luminosity would be 2 X 10^^ erg (corresponding to the length of the double arrows). This is in reasonable agreement with the 
observed values. 



intermediate region, mainly nucleon recoil in scattering 
reactions still allows for a certain amount of energy ex- 
change between the neutrinos and the medium, which 
reduces the mean energies of and Vt appreciably (Fig- 
ure[T5|). 

Our simulations show that this well-known effect on 
the spectra of f^/r has further consequences if the en- 
ergy exchange between the neutrinos and the background 
medium is included self-consistently. Neutrino absorp- 
tion and emission in the PNS surface region maintain 
a quasi-steady-state temperature stratification which is 
determined by the neutrino flux temperature and thus 
by the temperature of regions deeper in the neutron star 
(which is similar in models G15 and S15). In order to 
maintain this stratification, any additional energy input 
from i/^/t- at densities around 10^^ g cm~'^ must be com- 
pensated by a corresponding energy loss in and De- As 
this rather efficient conversion of i'^/t- "recoil energy" of 
the nucleons into Ve and De happens in a relatively hot 
layer, there is also a (minor) enhancement of the mean 
energies of Ve and Pe- We illustrate this "reshuffiing" of 
energy between the different flavors in Figure 1181 which 
shows that the luminosity decrease in I'fj./r due the scat- 
tering losses to nucleons and the enhancement of the i>e 
and De luminosity occur in precisely the same density 
region and correspond well in magnitude. 

5.3.4. Differences to Previous Pseudo-Newtonian 
Simulations 



Janka ( 200 9i) already inves tigated the 15Mq 
33 ()1995l ) in the framework 



Marek 

progenitor of .Wooslev fc Weaver 
of the effective potential approximation using the same 
ray-by-ray variable Eddington factor method that has 
been applied in our study. Given the close similarity of 



their model setup, the simulations presented here must 
also be interpreted against the background of their earlier 

r esults. 

iMarek &: Jankal (|2009[ ) considered four different sim- 
ulations of the 15Mq progenitor, namely two rotating 
models (128 angular zones) with different choices for 
the effective gravitational potential and different resolu- 
tion in energy space, and two highly-resolved runs (192 
angular zones) without rotation and with two different 
equations of state. Only one of these, the model includ- 
ing rotation and a gravitational potential that somewhat 
overe stimates strong-field effects (case R of Mare k et al.l 
120061) ■ developed an explosion roughly 550 ms after 
bounce. 

By combini n g ou r findings with those of 
IMarek fc Jallkal pOOl . we arrive at a rather co- 
herent picture. In both cases, the 15M0 progenitor 
appears to evolve very close to the explosion threshold, 
as small model variations were sufficient to bring about 
an explosion, while other simulations with slightly less 
favorable showed no signal of shock revival at least until 
> 400 ms after bounce. The qualitative simi l arities 
between our results and those of IMarek fc Janlal (|2009[ ) 
go much further: It is remarkable that in both cases the 
explosion is not linked to the transition of the Si/SiO 
interface through the shock but occurs at a much later 
stage when the SASI again becomes violent after an 
intermediate phase of relatively weak activity. In both 
cases, the most optimistic model is one with stronger 
neutrino e mission than for the "best" effective potential 
(case A) of IMarek et al.l (120061 ). which is either due to 
the better treatment of GR (for our model G15) or the 
choice of a stronger effective potential (case R) for the 
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explosion model LS-rot of iMarek fc Jankal ()2009[ fl. In 
the light of their earlier results, the beneficial role of GR 
corrections thu s seems all the m ore plausible. 

The work of IMarek fc Jankal (2009) also provides us 
with reference runs computed with a higher angular res- 
olution, and thus allow us to address an important limi- 
tation of the present first generation of GR neutrino hy- 
drodynamics models. In particular, they discuss a non- 
rotating effective potential model (LS-2D) corresponding 
exactly to model M15, but with a higher resolution of 192 
angular zones (as compared to 64). For the first 200 ms 
after bounce, the shock evolution, as well as the advec- 
tion and heating time-scal es agree remarkably wel l, while 
the high- resolution run of IMarek fc Jankal (|2009f) shows 
a somewhat more optimistic evolution afterwards. This 
suggests that our findings about the beneficial effects of 
GR will prove robust for two reasons: Differences be- 
tween the GR model G15 and the effective potential run 
M15 apparently assert themselves earlier than resolution 
effects and also appear to be more pronounced. More- 
ov er, the comparison of model M15 and model LS-2D 
of IMarek fc Jankal (j200 91 provides further evidence for 
the more optimistic evolution of hi gh-angular resolutio n 
runs in 2D that has been found bv iHanke et alJ (|2011l ). 
It therefore seems likely that high-resolution follow up 
studies will confirm our present results about explosions 
in GR. 

5.4. An Aside on Other Measures for the Criticality of 
the Accretion Flow 

Extending ear l ier work on the subjec t 
(iBurrows fc GoshV ImE; ' Yamasaki fc Yamadal 12001 . 
IPeicha jPrhompson (201^ recently re- investigated 
the problem of a stationary ID accretion flow onto a 
proto-neutron star in great detail, and attempted to 
derive firm criteria for the transition from the phase 
of quasi-stationary accretion to the explosion. They 
claimed to have found such a precise and robust criterion 
in the so-called antesonic condition for the maximum 
ratio of the local sound speed Cs and the escape velocity 
Vase hi^ the post-shock region, stating 



> 



3 
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(18) 



as requirement for an explosion. According to their re- 
sults, this antesonic condition is a more reliable measure 
for the criticality of the fiow than other criteria such as 
the ratio Tadv/Thoat used in SectionjSl or the requi rement 
of a g rowing mass a nd energy in the g ain region (jJankal 
120011 ). Very recently, iFernandezI (|201ir ) gave a very care- 
ful assessment of the Burrows-Goshy-limit and its rela- 
tion to the conditions where a steady-state accretion fiow 
in ID meets the threshold for a runaway instability. 

Considering the ongoing debate about explosion cri- 
teria, it is worthwhile to study the applicability of the 
new antesonic condition to dynamical multi-dimensional 
simulations of the accretion phase, for which the assump- 

^ The explosion model LS-rot of Marck & Janka ( 2009) also in- 
cludes rotation, which, as these authors found, adversely affects the 
heating conditions. This mitigates the effect of the overly strong 
gravitational potential and makes their model explode at a sim- 
ilarly late stage as model G15, although with an extra delay of 
> 100 ms. 
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Figure 19. Time evolution of the maximum value of the ratio of 
the sound speed Cs and the (local) escape velocity Desc for models 
Gil, G15, M15, S15, and N15, smoothed over 10 ms. Two different 
methods for evaluating max{c^ / v^^^) are used: We either compute 
spherical averages of before taking the maximum (case A) or 
take the global maximum (case B). Only the region outside the 
PNS (p < 10^^ g cm""^) is considered in either case. The escape 
velocity is computed directly from the gravitational potential or, 
in the GR case, the lapse function. 



tions of IPeicha fc ThompsonI ()2012D (e.g. spherical sym- 
metry, stationarity) do not necessarily constitute a valid 
approximation. In order to gauge the usefulness of the 
new criterion in this context, we need to address two 
specific questions: First, how sharply does the antesonic 
condition resolve differences between optimistic and pes- 
simistic models in the accretion phase, and, second, can it 
accurately capture the onset of the runaway and the de- 
velopment of an explosion? To this end, we plot the time 
evolution of max{c'^ / v^^^) for our 2D models in Figure 1191 
using two different evaluation methods, viz. averaging 
Cs over solid angle before taking the maximum (Case A), 
and directly computing the maximum of the ratio cj. /f osc 
in the post-shock region (Case B). As the global maxi- 
mum is very sensitive to fluctuations and not necessarily 
indicative for the overall conditions in the gain region. 
Cas e A seems more in keeping with the ID analysis of 
'Pejcha fc Thompson (2012). 

Interestingly, we do not observe any clear differences in 
max(c^/Wpg(.) between the 2D models in Case A prior to 
the development of an explosion. After an initial growth 
from 0.1 to 0.17, max(c^/D^g^) settles to a stable average 
value of about 0.17 in all our simulations of the 15Mq 
progenitor. Furthermore, the evolution of max(c^/wpj,j,) 
in ID is also virtually indistinguishable from 2D in the 
case of the non-exploding model S15 (not shown in Fig- 
ure [in] as the curves would almost overlap) If we 
compute the absolut e maxi mum of the velocity ratio of 
IPeicha fc ThompsonI ()2012D between the proto-neutron 
star and the shock (case B) , we see an increase to values 
between 0.2 and 0.35 once violent convective and SASI 
activity starts (as a result of the higher entropy in the 
convective plumes and SASI lobes). However, a clear 
hierarchy between the different 15Mq models is again 
absent; model M15, which appears to be the most pes- 

* High values of raax(c1 / v'^^^) > 0.15 are found in all our ID 
models. Such high values naturally follow from the jump con- 
ditions at the shock for the typical infall velocities in dynami- 
cal simu lations, which are considerably higher than assumed by 
IPeicha Thompson (2012). It is unclear to what extent the choice 
of the infall velocity or the omission of GR affects the findings of 
IPeicha fc 'I'hompsonI |[20T2'| . 
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simistic case according to the time-scale ratio Tadv/'T'heat, 
even reaches somewhat higher values than model S15 on 
average. 

On the other hand, the onset of the explosion is in- 
deed correlated with an increase in ra&yi(c1 / v^^^ beyond 
0.2. In Case A this increase can be delayed signifi- 
cantly for a strongly asymmetric explosion (as occurs in 
run G15). The absolute maximum of c^/Wesc appears 
to be a more robust indicator, but the critical value 
would have to b e quite different from that advocated by 
iPejcha k, Thom pson (2012), namely, 

max {cl/vi:) > 0.35, (19) 

but the evidence from two explosion models and two 
non-exploding models of one of the progenitors is far 
from conclusive. We find no compelling link between 
this (incidental?) observa tion and the arguments of 
iPeicha fc ThompsonI (|20T2f ). 

Thus, unlike the time-scale ratio Tadv/Theat, the an- 
tesonic condition does not appear to be a very use- 
ful measure for distinguishing optimistic and pessimistic 
models during the accretion phase, and is of limited 
use for detecting the runaway situation that ultimately 
leads to the explosion. There is a number of possible 
reasons for this app a rent c onflict with the findings of 
IPeicha fc ThompsonI (|2012f) . Different from the time- 
scale ratio Tadv/Theat, the antesonic condition does not 
make any explicit reference to the heating conditions in 
the gain layer; it depends only on the pressure, density 
and gravitational potential at a specific point in the ac- 
cretion flow (corresponding more or less to the gain ra- 
dius), which may be crucial for the existence of a sta- 
tionary accretion flow in spherical symmetry. What hap- 
pens to such a local quantity beyond the critical point 
in the dynamical flow of a multi-dimensional situation is 
unclear, however, whereas Tadv/Theat contains global in- 
formation directly related to the competing influence of 
accretion and neutrino heating. 

Moreover, convection and the SASI alter the structure 
of the post-shock accretion flow, the heating conditions, 
and the propagation of the shock considerably in a multi- 
dimensional setup. The stratification of the gain layer 
is not only modified by turbulent mixing; the deviation 
of the velocity field from that of a stationary spherical 
accretion fiow and the deformation of the shock also be- 
come large enough to invalidate many of the assump- 
tions the one-dimensional model of IPeicha fc ThompsonI 
(|2012D - such as the time-independent, spherically sym- 
metric Euler equations (neglecting, e.g., the turbulent 
pressure) and the simple form of the jump conditions 
at the shock (neglecting the SASI). More extreme cases 
with simultaneous accretion and shock expansion in a 
dipolar explosion (such as model G15) can hardly be 
accommodated in a ID setup at all. These complica- 
tio ns certainly have the potent ial to change the findings 
of IPeicha fc ThompsonI ()2012D both quantitatively and 
qualitatively, and to make meLx{c1 / v"^^^) an inferior run- 
away criterion for multi-dimensional supernova simula- 
tions. As demonstrated in the preceding sections, the 
time-scale ratio Tadv/Theat is less affected by such com- 
plications and therefore remains, in our opinion, a better 
diagnostic quantity for the evolution towards the explo- 
sion. 



6. SUMMARY AND CONCLUSIONS 

In this paper, we have presented the first multi- 
dimensional (2D) general rclativistic simulations of su- 
pernova explosions with a sophisticated neutrino trans- 
port treatment. U sing the recently in troduced Vertex- 
CoCoNuT code ()Miiller et al.ll2010l) . which is based on 
a rclativistic generalization of the va riable Eddington 
factor method for neutrino transport (jRampp fc Jankal 
|2002| ) combined with the "ray-by -ray-plus" approach 
for multi-dimension al problems ()Buras et al.l l2006bl : 
iBruenn et all l2006l) . we calculated the evolution of 
two pro genitors with ll.2Mr:~, (IWooslev et al.ll2002D and 
ISMq (jWooslev fc Weaver] [1995I ). For the 15Mg case, 
we also performed complementary simulations with a 
pseudo-Newtonian and a purely Newtonian treatment of 
gravity, and with a simplified set of neutrino interaction 
rates. We have conducted a comprehensive analysis of 
the shock propagation, the explosion dynamics, and the 
heating conditions for these models. With regard to the 
major goals formulated in Section [1] our results can be 
summarized as follows: 

1. Feasibility of multi-dimensional rclativistic super- 
nova simulations: For both the 11.2M0 and the 
15Mq progenitor, we obtain explosions in 2D and 
have been able to extend the simulations several 
hundreds of milliseconds into the post-bounce and 
the explosion phase, thus proving the rclativistic 
Vertex-CoCoNuT code to be as robust and sta- 
ble as its Newtonian counterpart Prometheus- 
Vertex. 

2. Qualitative verification of pseudo-Newtonian re- 
sults: For the 11.2 M(^ progenit or, we con- 
firm the explos i ons o f iBuras et al.l (|2006a) and 
iMarek fc Ja"Skal (|2009D . and for the 15Mq pro- 
genitor, the evolution is very similar to differ- 
ent models computed with the pseudo-Newtonian 
Prometheus- Vertex code for the first ~ 
400 ms (i.e. before the onset of the explosion 
in our GR model). As Vertex-CoCoNuT and 
Prometheus- Vertex employ completely inde- 
pendent hydro solvers and somewhat different rou- 
tines for the moment equations, our results can 
be viewed as an important step towards code ver- 
ification in the sense that different codes based 
on a similar physical model have been shown to 
produce similar results. However, the different 
physics leads to important quantitative differences 
(see next item). 

3. Role of general relativistic effects in the su- 
pernova problem: For t he 15Mq progenitor of 
IWooslev fc Weaver! (|1995| ). we performed a com- 
parison of a relativistic, a purely Newtonian, and 
a pseudo-Newtonian simulation in order to de- 
termine the importance of GR effects during the 
post-bounce phase. We find significant differences 
on the order of several tens of percent in some 
explosion-relevant quantities with a tendency to- 
wards more optimistic heating conditions in the GR 
case. This diagnosis is confirmed by the fact that 
we observe an incipient explosion in the GR run 
at ~ 400 ms, whereas the Newtonian and pseudo- 
Newtonian models fail to explode, at least until 
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~ 500 ms after bounce. Different thermodynamic 
conditions in the neutrinospheric region have been 
identified as the most likely cause for the more 
optimistic evolution in GR. It is noteworthy that 
the pseudo-Newtonian approach - despite a much 
better quantitative agreement with the GR model 
than for the Newtonian run - still leads to a dif- 
ferent outcome than the GR case. Interestingly, 
our GR result yields a similar (successful) lat e ex- 
plosion as obtained bv lMarek fc JanO (|2009l ) or a 
pseudo-Newtonian 15Mq model including rotation 
and a gravitational potential that has been found 
to somewhat overestimate GR corrections in spher- 
ical symmetry. 

Thus, contrary to claims in the literature 
(jNordhaus et al.l[2010( ) , general relativity may class 
as a "tens of percent" effect with a similar impact 
on the dynamics as dim ensional itv ([Nordhaus et al.l 
I2010I: iHanke et al.llM l: Takiw aki et al.ll20lH) and 
the equation of state (Marek et al., in preparation). 
GR thus emerges as a key ingredient for accurate 
supernova simulations. 

4. Influence of the neutrino physics input: Based on 
a 15Mq model computed with a simplified set of 
opacities, we also illustrated the sensitivity of the 
neutrino heating conditions to the neutrino inter- 
action rates. Specifically, we found that the energy 
transfer from v^/t to the background medium in 
the neutrinospheric region by nucleon recoil has 
a beneficial effect, since it effectively provides a 
means of "converting" a few 10^^ erg s~^ from 
v^/t vcAo Ve and Pe, which can then actively con- 
tribute to the heating in the gain layer. Without 
the improved opacities, we observe at least a con- 
siderable delay of the explosion in the I^Mq case 
(to > 450 ms after bounce). Our findings further 
substantiate claims (jBruenn et al.| [2b09) about the 
importance of state-of-the-art neutrino physics in- 
put for supernova simulations, although the effect 
appears to be smaller than that of GR. 

On the whole, the inclusion of general relativity in su- 
pernova models has turned out to be much more than 
a marginal improvement, even compared to the effec- 
tive potential approach which has hitherto been used by 
some groups (Rampp fc Janka 2002: M a rek et al. 2006; 
IBruenn et a"ni2009l: IScheidegger et al.ll2008D . The purely 
Newtonian approach is definitely ruled out as a viable ba- 
sis for quantitatively accurate models of core-collapse su- 
pernovae. GR may prove more important than expected 
for accurately capturing the physics of neutrino-driven 
explosions if the 15Afo progenitor is any indication. In 
a subsequent publication, we will demonstrate that GR 
also has a large impact on the observational signatures 
such as neutrinos and gravitational waves from super- 
nova cores. 

This adds another interesting facet to the supernova 
problem and probably indicates that the explosion mech- 
anism does not hinge on a single dominating factor af- 
ter all. General r elativity, neutrino microphysics, the 
equation of sta te (iMarek fc Janka 200 !^, and dimen- 
sionality effects (N ordhaus et al.1120101: IHanke et al.ll2011l: 
pTakiwaki et al..,20lll ) mav eauallv contribute to fill the 



missing parts of the current picture. 

Such a situation naturally opens several avenues for 
future research efforts. Above all, the potentially ben- 
eficial effects of GR need to be confirmed by high- 
resolution studies (although improved angular resolution 
will, in all prob ability, only lead t o more robust explo- 
sions in 2D, cp. IHanke et al.ll20lTl ). In the long run, a 
comprehensive and self-consistent approach to supernova 
modeling, covering the essential aspects of the problem 
(such as neutrino transport, general relativity, 3D hy- 
drodynamics) with the help of highly accurate numeri- 
cal algorithms, will be indispensable for a firm quantita- 
tive understanding of core-collapse supernova explosions. 
Vertex-CoCoNuT provides a possible platform for fur- 
ther developments in that direction, e.g. the extension of 
the code to 3D, or the inclusion of more accurate multi- 
angle transport. On the other hand, a more thorough 
understanding of the individual components of the su- 
pernova problem (interplay of SASI and convection, crit- 
ical explosion conditions, nucleosynthesis, neutron star 
kicks and spins, etc.) is equally important, and this is 
no less true for the effects of general relativity discussed 
in this paper. Having confined ourselves to two progen- 
itors with 11.2A/0 and 15M0, which produce neutron 
stars of moderate compactness, we believe that one of 
the essential tasks will be that of probing deeper into 
the strong-field regime, where general relativity can be 
expected to play an even greater role. In our view, the 
case of more massive progenitors, possibly on the verge 
to black hole formation, therefore merits particular at- 
tention in the future. 
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